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ABSTRACT 

Graph-based spectral denoising is a low-pass filtering using 
the eigendecomposition of the graph Laplacian matrix of a 
noisy signal. Polynomial filtering avoids costly computa¬ 
tion of the eigendecomposition by projections onto suitable 
Krylov subspaces. Polynomial filters can be based, e.g., 
on the bilateral and guided filters. We propose constmct- 
ing accelerated polynomial filters by running flexible Krylov 
subspace based linear and eigenvalue solvers such as the 
Block Locally Optimal Preconditioned Conjugate Gradient 
(LOBPCG) method. 

Index Terms — Image denoising, spectral polynomial fil¬ 
ter, graph Laplacian, Krylov subspace method 

1. INTRODUCTION 

In this note, we deal with noise removal from a given noisy 
signal, which is a basic problem in signal processing, with ap¬ 
plications, e.g., in image denoising m. Apart from the trivial 
application of removing noise prior to presenting the image to 
a human observer, pre-smoothing an image and noise removal 
may help to improve performance of many image-processing 
algorithms, such as compression, enhancement, segmentation 
etc. A noise removal operation is often referred to as a filter. 

Modern denoising algorithms endeavor to preserve the 
image details while removing the noise. One of the most pop¬ 
ular denoising filters is the bilateral filter (BF), which smooths 
images while preserving edges, by taking the weighted av¬ 
erage of the nearby pixels. The weights depend on both the 
spatial distance and photometric similarity, thus providing 
local adaptivity to the input image. Bilateral filtering was 
introduced in 013] SI as an intuitive tool without theoret¬ 
ical justification. Since then, connections between the BF 
and other well-known filtering techniques such as anisotropic 
diffusion, weighted least squares, Bayesian methods, kernel 
regression and non-local means have been explored, see, e.g., 
survey 0. 

A convenient way to represent images is by graphs 0, 
especially when considering images over irregular grids and 
treating non-local interactions between pixels. A single appli¬ 
cation of the bilateral filter to an image may be interpreted as 


a vertex domain transform on a graph with pixels as vertices, 
intensity values of the pixels as the graph signal and filter co¬ 
efficients as link weights that capture the similarity between 
the vertices. The BF transform is a nonlinear anisotropic dif¬ 
fusion, cf. QElIll, determined by entries of the graph Lapla¬ 
cian matrix, which are related to the BF weights. In the lin¬ 
ear case, the solutions of time-dependent anisotropic diffusion 
problems are represented by operator exponentials or semi¬ 
groups, which can be approximated by operator polynomials. 

Eigenvalues and eigenvectors of the graph Laplacian ma¬ 
trix allow us to apply the Fourier analysis to the graph signals 
or images as in 0 and perform frequency selective filtering 
operations on graphs, similar to those in traditional signal pro¬ 
cessing. Expensive computation of eigenvectors is avoided 
in polynomial spectral filters, which are fully implemented in 
the vertex domain. Eor example, the spectral filters in Gojini 
are based on optimal polynomials constructed by means of the 
Chebyshev approximation and conjugate gradient algorithm. 

A recent newcomer in the field of filtering is the guided 
filter (GE) proposed in f[2[ [131 and included into the image 
processing toolbox of MATLAB. According to our limited 
experience, the GE filter is significantly faster than the BE 
filter. The authors of ca additionally advocate that GE is 
gradient preserving and avoids the gradient reversal problems 
in contrast to BE, which is not gradient preserving. 

Iterative application of smoothing filters like BE and GE 
can be interpreted as matrix power transforms, in general 
case, nonlinear, or, equivalently as explicit integration in 
time of the corresponding anisotropic diffusion equation. In 
the present paper, we continue the approach of nnidi] and 
propose acceleration of the iterative smoothing filters based 
on the polynomial approximations implicitly constructed in 
the conjugate gradient (CG) algorithm and in the LOBPCG, 
which is a leading eigensolver for large symmetric matri¬ 
ces |[l4l . LOBPCG has been successfully applied to image 
segmentation in ca. 

The remainder of the paper is organized as follows. Sec¬ 
tion 2 gives main formulas of the original bilateral filter in 
the graph-based notation. Section 3 provides a similar de¬ 
scription of the guided filter from na. Section 4 introduces 
a basic Eourier calculus on graphs and applies it to a sim¬ 
plified frequency analysis of iterations with the graph Lapla- 
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cian matrix. Section 5 describes the CG acceleration of the 
smoothing filters. Section 6 gives a brief description of the 
LOBPCG method and suggests it as an acceleration for BF 
and GF. Our numerical experiments in Section 7 demonstrate 
improved performance of our accelerated hlters in the one¬ 
dimensional case mostly for the BF hlter. Similar improve¬ 
ments are expected when accelerating the GF hlter. 


second way is also faster to compute since the BF weights are 
computed only once in the beginning. 

The fastest implementations of BF have the arithmetical 
complexity 0{N), where N is the number of pixels or voxels 
in the input image uainiiiii. 

3. GUIDED FILTER AS A VECTOR TRANSFORM 


2. BILATERAL FILTER AS A VECTOR 
TRANSFORM 


The bilateral hlter transforms an input image x into the output 
image y by the weighted average of the pixels of x: 


y^ = 


E 


Wij 


Lj. 


( 1 ) 


Let Pi denote the geometrical position of a pixel i. Then 


Wij = exp 


\\Pi-P3\\^ \ 

2^? J 


exp 



2a2 J ’ 


( 2 ) 


where CTs and ar are the hlter parameters JJ), \\pi — pj\\ is 
a spatial distance between pixels i and j. For color images, 
the photometric distance \xi — Xj\ can be computed in the 
CIE-Lab color space as suggested in a. 

The BF weights Wij determine an undirected graph G = 
(V, £), where the vertices V — {1,2,..., N} are the pixels 
of the input image and the edges S = {{i,j)} connect the 
“neighboring” pixels i and j. The adjacency matrix W of the 
graph G is symmetric and has the entries Wij > 0. Let D 
be the diagonal matrix with the nonnegative diagonal entries 
di — In the matrix notation, the bilateral hlter op¬ 

eration (fU is the vector transform with W depending on a 
guidance image (we use a; as a guidance image in (|^) 

y = D-^Wx = x-D-^Lx, (3) 


Algorithm I Guided Filter 

Input: X, g, p, e 

Output: y 

UlCCLTig = fmean(^9^P) 

TTlCUnx — fmeani^T P') 

CorVg = fmeanig- * g, p) 

COrVgx = fmeanig- * X, p) 

vaxg = corVg — meaug. * meaug 

covgx = corvgx — meaUg. * mearix 

a = covgx-livavg -f e) 
b = meaUx — a. * meaUg 
meaUa = fmean{a,p) 
mearib = fmeanib, p) 

y = meana- * g + mearii, 

Algorithm 1 is a pseudo-code of the guided hlter from 
iia, where fmeani', p) denotes a mean hlter with the win¬ 
dow width p. The constant e determines the smoothness de¬ 
gree of the hlter: the larger e the larger a smoothing effect. 
The arithmetical operations .* and ./ are the componentwise 
multiplication and division. 

The input image is x, the output image is y. The guidance 
image g has the same size as x. If g coincides with x, then 
the guided hlter is called self-guided. The arithmetical com¬ 
plexity of GF equals 0{N), where N is the number of pixels 
in X. 

The guided hlter operation of Algorithm 1 is the transform 
y = W{g)x, (5) 


where 

L = D-W. (4) 

is the Laplacian matrix for the graph G. Gershgorin’s theo¬ 
rem from matrix analysis guarantees a smoothing effect of the 
bilateral hlter, when max^ di < 2. 

The bilateral hltering can be used iteratively. There are 
two ways to iterate BF: (1) by changing the weights Wij at 
each iteration using the result of the previous iteration as a 
guidance image, or (2) by using the hxed weights at each iter¬ 
ation as calculated from the initial image as a guidance image 
for all iterations. The hrst alternative results in a nonlinear hl¬ 
ter, where the BF graph changes at every iteration and we can 
only provide separate spectral interpretation for each stage. 
The second alternative generates a linear hlter, the Laplacian 
matrix is hxed for all iterations, and it is possible to provide a 
spectral interpretation of the whole cascaded operation. The 


where the implicitly constructed symmetric matrix W(g) has 
the following entries, see ini: 


W^jig) 



E 
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L, idi- Pk)ig3 
\ <xl + e 



(6) 

Here ujk is the neighborhood around pixel k of width p, where 
the mean hlter fmeani‘,p) is applied, and |a;| denotes the 
number of pixels in ojk, the same for all k. The values pk 
and af. are the mean and variance of the image gmu3k- 

The graph Laplacian matrix is obtained from the matrix 
W in the standard way. According to ifT^ . the values di = 
Wij equal 1. Therefore, the diagonal matrix D equals the 
identity matrix I, and the graph Laplacian matrix is the sym¬ 
metric nonnegative dehnite matrix L = I — W. 

The guided hltering can also be used iteratively. Similar 
to the BL hlter, the guided hlter is applied iteratively: (1) 








either using the result of the previous iteration as a guidance 
image, or (2) by using the initial image as a guidance image 
for all iterations. The hrst alternative results in a nonlinear 
GF hlter. The second alternative generates a linear GF hlter, 
and the Laplacian matrix remains hxed at each iteration. 

4. SPECTRAL INTERPRETATION OF THE 
LOW-PASS FILTERS 

The graph eigenstructure is the eigenvalues and eigenvectors 
of the graph Laplacian matrix L, which is symmetric and non¬ 
negative dehnite. In certain situations, the normalized Lapla¬ 
cian matrix (diagL)“^/^L(diag may be more suitable 

than L. The spectral factorization of L is the matrix decom¬ 
position 

L = UAU'^, (7) 

where the diagonal elements Ai of the diagonal matrix A and 
columns Ui of the orthogonal matrix U = [rti,. .., u„] are, 
respectively, the eigenvalues of L and corresponding eigen¬ 
vectors of the unit 2-norm. 

Similar to the classical Fourier transform, the eigenvec¬ 
tors and eigenvalues of the graph Laplacian matrix L provide 
the oscillatory structure of graph signals. The eigenvalues 
Ai < • • • < Ai < • • • < A„ can be treated as graph fre¬ 
quencies. The corresponding eigenvectors Ui of the Lapla¬ 
cian matrix L are generalized eigenmodes and demonstrate 
increasing oscillatory behavior as the magnitude of the graph 
frequency increases. The Graph Fourier Transform (GFT) of 
an image x is dehned by the matrix transform x = U"^x, the 
inverse GFT is the transform x = Ux. 

Let us consider the BF vector transform Q. In numeri¬ 
cal analysis, the linear transformations are called 

the power iterations with the amplihcation matrix D~^W or 
simple iterations for the equation Lx — 0 with the precon¬ 
ditioner D. Application of the transform Q preserves the 
low frequency components of x and attenuates the high fre¬ 
quency components, cf. nniiii]. It is also well-known that 
the Krylov subspaces well approximate the eigenvectors cor¬ 
responding to the extreme eigenvalues. Thus the projections 
onto the suitable Krylov subspaces would be an appropriate 
choice for high- and low-pass biters. The Krylov subspace 
methods are efficient owing to their low cost, reasonably good 
convergence and simple implementation without painful pa¬ 
rameter tuning. The convergence can be accelerated by the 
aid of good preconditioners. 

5. CG ACCELERATION 

Since the graph Laplacian matrix L is symmetric and non¬ 
negative dehnite, the hrst candidate to replace the simple iter¬ 
ation Xk+i = Xk — D~^Lxk is the preconditioned conjugate 
gradient (CG) method for the system of homogeneous linear 
equations Lx = 0 with the preconditioner matrix D. To avoid 


over-smoothing, only few iterations of the preconditioned CG 
method must be performed. A large number of iterations for 
Lx — 0 produces piecewise constant images and is better 
applicable to the image segmentation problems. Good refer¬ 
ences for the Krylov subspace methods are the books 

Algorithm 2 Preconditioned CG(fcmax) for Lx = 0 

Input: L, Xq, fcmax, preconditioner D 

Output: X 

X = Xq 

r = —Lx 

for /c = 1, . . . , fcmax do 

s = D~^r 

if fc = 1 then 

p = s 

else 

j3 = (s^(r - roid))/{sMSoid) 

p = s + (3p 

endif 

q = Lp 

a = {s^r)/(p'^q) 

X = X -\- ap 

Told = r 

Sold = s 

r = r — aq 

endfor 

Algorithm 2 is a slightly modibed standard precondi¬ 
tioned CG algorithm formally applied to the system of linear 
equations Lx = 0. It contains the formula for /3 that is 
different from that in the MATLAB implementation of the 
preconditioned CG. Such a formula converts the algorithm 
into a bexible variant, which possesses better convergence 
properties, when the input matrices L and D may change; 
see, e.g., 112111221 . 

The CG iterations are the polynomial biters, i.e., repre¬ 
sented in the form Xk = Pk{L)xQ, where the coefficients of 
the polynomial pk (A) of degree k depend on the matrix L and 
initial vector xq. In the preconditioned case, the biter has the 
form Xk = pk{D~^L)xQ, and the coefficients of the polyno¬ 
mial p/j(A) depend on D~^L and xq. 

6. LOBPCG ACCELERATION 

A more powerful alternative to the Krylov subspace solvers 
for linear systems as the polynomial spectral biters are the 
eigensolvers with preconditioning. Since the graph Laplacian 
matrix L is symmetric, we propose the symmetric eigensolver 
LOBPCG with preconditioning for construction of low-pass 
biters. LOBPCG is a very efficient eigensolver and gives 
a fast solution to the spectral image segmentation problem, 
see ca. To avoid over-smoothing during noise removal, 
LOBPCG must perform only few iterations. 

Algorithm 3 below is a standard non-blocked version of 
the preconditioned LOBPCG algorithm, which smooths the 




input signal xq with respect to the eigenmodes of the Lapla- 
cian matrix L. If necessary, Algorithm 3 can be modihed to 
obey the constraint that the vectors Xk must be orthogonal to 
the vector e with all components equal to 1. The vector e is 
the eigenmode corresponding to the zero eigenvalue of L. 

Algorithm 3 LOBPCG method as a low-pass hlter 

Input: L, D, xq, and a preconditioner T 

Output: Xkmax 
Po = 0 

for fc = 0, . . . , /Cmax - 1 do 
Xk = {xl Lxk) / {xl Dxk) 

r = Lxk - XkDxk 
Wk = Tr 

use the Rayleigh-Ritz method for the pencil L — XD 
on the trial subspace span{wk, Xk,Pk} 

Xk-\-l — tJJk TkXk T ^kPk 

(the Ritz vector for the minimum Ritz value) 

Pk+l =Wk+ 'JkPk 
endfor 

7. NUMERICAL EXPERIMENTS WITH THE 
KRYLOV SUBSPACE-BASED POLYNOMIAL 
ACCELERATION OF THE FILTERS 

Our experiments have been done in MATLAB. We have com¬ 
pared performance of the bilateral hlter versus the CG ac¬ 
celerated BF (BF-CG) hlter versus the self-guided GF hlter 
for several one-dimensional signals. The clean ID signal is 
piecewise linear with steep ascents and descents. The noisy 
signal is obtained from the clean signal by adding the Gaus¬ 
sian noise 0.1 • randn{N, 1), where we test N = 500 in Fig¬ 
ure [T] and N = 1000 in Figure where N is number of 
samples in the signal discretized on a uniform grid. 

Figure [T] and Figure [^display our test results comparing 
500 iterations of the self-guided BF with the window width 
equal to 3, 20 iterations of the self-guided GF with the MAT- 
LAB default window width 5, and, hnally, 20 iterations of the 
CG accelerated BF, using the noiseless signal as a guidance 
and the diagonal matrix D as the preconditioner. 

In Figures and 1^ we observe high quality denoising by 
all tested hlters, both in terms of the average ht of the noise¬ 
less signal and of preserving the signal details. We notice that 
GF has a tendency to round sharp corners compared to BF and 
BF-CG, which is not surprising, since GF performs extensive 
smoothing on relatively wider neighborhoods. Interestingly, 
GF error is strongly discontinuous on the signal edges, while 
BF and BF-CG provide more conservative approximations. 
Increasing the number of samples apparently makes the BF 
and BF-CG errors smoother as can be seen comparing Fig¬ 
ure [T] and Figure]^ 

Miraculously all three hlters give nearly the same error 
with the reported parameters, which requires further investi¬ 
gation and explanation. 




Fig. I. Filters in ID case; 500 points 
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Fig. 2. Filters in ID case; 1000 points 


Our second series of tests deals with image denoising. We 
use the 5-point stencil to generate the lattice for the graph 
Laplacian matrix and the edge weights are determined by the 
BF hlter with = 0.1. 

Twenty iterations are performed using the CG and LOBPCG 
accelerated hlters with and without the constraint. The quality 
of denoising is similar improving PSNR from approximately 
20 for the noisy image to approximately 21 for the hltered 
images. The edges are well preserved but one can notice 
forming salt and pepper noise, which gets more extensive 
if the number of iterations is increased. Figure compares 
the BF hlter with the window half-width 5 and bilateral hlter 
with the standard deviation a = 0.1. 







































Fig. 3. Noisy versus clear image: gaussian noise, mean = 0, 
variance = 0.01, PSNR = 20.1 



Fig. 4. BF filter versus CG 


Conclusions 

We propose accelerating iterative smoothing filters such as the 
bilateral and guided filters, by using the polynomials, gener¬ 
ated by the conjugate gradient-type iterative solvers for lin¬ 
ear systems and eigenvalue problems. This results in efficient 
parameter-free polynomial filters in the Krylov subspaces that 
approximate the spectral filters corresponding to the graph 
Laplacians, which are constructed using the guidance signals. 
Our numerical tests for one-dimensional signals demonstrate 
the typical behavior of the proposed accelerated filters and 
explain the motivations behind the construction. The tests 



Fig. 5. LOBPCG with versus without the constraint 


Fig. 6. LOBPCG with the constraint versus CG 

showing image denoising reveal that the proposed filters are 
competitive. Our future work will concern accelerating the 
guided filters, testing the influence of non-linearities on the 
filter behavior, and incorporating efficient preconditioners. 
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